/*
 * Cerenkov_Counter3.cc
 *
 * Created: October 2012
 * Authors: Jessie McNew, David Blyth
 */

#include <cmath> //
#include <sstream>
using namespace std;

#include <G4Box.hh>
#include <G4Color.hh>
#include <G4LogicalBorderSurface.hh>
#include <G4LogicalSkinSurface.hh>
#include <G4OpticalSurface.hh>
#include <G4PVPlacement.hh>
#include <G4SDManager.hh>
#include <G4VisAttributes.hh>

#include "RICHPMTDetector.hh"
#include "Cerenkov_Counter3.hh"

Cerenkov_Counter3::Cerenkov_Counter3() {
	// Build materials

	nistMan = G4NistManager::Instance();

	Al = nistMan->FindOrBuildMaterial("G4_Al");

	SiO2 = nistMan->FindOrBuildMaterial("G4_SILICON_DIOXIDE");

	double airEnergy[2] = { 2.08 * eV, 5.90 * eV };
	double airRefIndex[2] = { 1.0003, 1.0003 };
	double airAbsorp[2] = { 1000 * m, 1000 * m };

	air = nistMan->FindOrBuildMaterial("G4_AIR");
	G4MaterialPropertiesTable* airMPT = new G4MaterialPropertiesTable();
	airMPT->AddProperty("RINDEX", airEnergy, airRefIndex, 2);
	airMPT->AddProperty("ABSLENGTH", airEnergy, airAbsorp, 2);
	air->SetMaterialPropertiesTable(airMPT);

	// Other initializations
	airPhys = NULL;
}

Cerenkov_Counter3::~Cerenkov_Counter3() {
	;
}

void Cerenkov_Counter3::Build(G4VPhysicalVolume* worldPhys) {
	// *** Array element positions correspond to unique tiles. ***

	const int M = 21; // array length
	const int N = 3; // Number of tiles being used
	double aerogelLength[N] = { 10 * cm, 10 * cm, 10 * cm };
	double aerogelWidth[N] = { 10 * cm, 10 * cm, 10 * cm };
	double aerogelDepth[N] = { 1.6 * cm, 1.6 * cm, 1.6 * cm };
	double aerogelZ[N] = { 24.3 * cm, 22.2 * cm, 20.6 * cm }; //todo: replace Z-coordinates w/ G4Vectors

	double aerogelGreenIndex[N] = { 1.021, 1.025, 1.029 };

	double PMTGridWidth = 10 * cm;
	double PMTGridLength = 10 * cm;
	double PMTGridDepth = 1 * cm;

	double PMTGridOpticalEnergy[M] = { 6.1992, 5.5104, 4.9594, 4.5085, 4.1328,
			3.8149, 3.5424, 3.3062, 3.0996, 2.9173, 2.7552, 2.6102, 2.4797,
			2.3616, 2.2543, 2.1562, 2.0664, 1.9837, 1.9074, 1.8368, 1.7712 };
	for (int i = 0; i < M; i++)
		PMTGridOpticalEnergy[i] *= eV;
	double PMTGridOpticalRefl[M] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
			0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
	double PMTGridOpticalEff[M] = { 0, 0, .01, .116, .240, .318, .332, .340,
			.332, .318, .263, .222, .163, .100, .060, .040, .024, .016, .004,
			.0, .0 };
	for (int i = 0; i < M; i++)
		PMTGridOpticalEff[i] *= 1;
//	double PMTGridOpticalEff[M] = { .090, .125, .163, .218, .228, .318 * .77,
//			.332 * .77, .340 * .77, .332 * .77, .318 * .77, .263 * .77, .222
//					* .77, .163 * .77, .100 * .77, .060 * .77, .040 * .77, .024
//					* .77, .016 * .77, .004 * .77, .0, .0 };

// Build optical surfaces

	G4OpticalSurface* PMTGridOptical = new G4OpticalSurface(
			"PMTGrid_photocathodeOptical", glisur, polished, dielectric_metal);
	G4MaterialPropertiesTable* PMTGridOpticalMPT =
			new G4MaterialPropertiesTable();
	PMTGridOpticalMPT->AddProperty("REFLECTIVITY", PMTGridOpticalEnergy,
			PMTGridOpticalRefl, M);
	PMTGridOpticalMPT->AddProperty("EFFICIENCY", PMTGridOpticalEnergy,
			PMTGridOpticalEff, M);
	PMTGridOptical->SetMaterialPropertiesTable(PMTGridOpticalMPT);

// Build volumes

	RICHPMTDetector* cSensitiveDetector = new RICHPMTDetector("RICH/PMT");
	G4SDManager* SDMan = G4SDManager::GetSDMpointer();
	SDMan->AddNewDetector(cSensitiveDetector);

	G4Box* airSolid = new G4Box("Air_Solid", 50 * cm, 50 * cm, 50 * cm);
	G4LogicalVolume* airLog = new G4LogicalVolume(airSolid, air, "Air_Log", 0,
			cSensitiveDetector);
	airLog->SetVisAttributes(G4VisAttributes(false, G4Color::Yellow()));
	airPhys = new G4PVPlacement(0, G4ThreeVector(0, 0, 0), "Air_Phys", airLog,
			worldPhys, false, 0);
	(void) airPhys;

	G4Box* PMTGridSolid = new G4Box("PMTGrid_Solid", 0.5 * PMTGridWidth,
			0.5 * PMTGridLength, 0.5 * PMTGridDepth);
	G4LogicalVolume* PMTGridLog = new G4LogicalVolume(PMTGridSolid, Al,
			"PMTGrid_Log");
	PMTGridLog->SetVisAttributes(G4VisAttributes(true, G4Color::Black()));
	G4LogicalSkinSurface* PMTGridSurface = new G4LogicalSkinSurface(
			"PMTGrid_photocathodeSurface", PMTGridLog, PMTGridOptical);
	(void) PMTGridSurface;
	G4VPhysicalVolume* PMTGridPhysical = new G4PVPlacement(0,
			G4ThreeVector(0 * cm, 0 * cm, -PMTGridDepth / 2.), "PMTGrid_Phys",
			PMTGridLog, airPhys, false, 0);
	(void) PMTGridPhysical;

	for (int i = 0; i < N; i++) {
		stringstream ss;
		ss << "Aerogel" << i << "_Log";
		G4LogicalVolume* aerogelLog = buildAerogelTile(ss.str().c_str(),
				aerogelLength[i], aerogelWidth[i], aerogelDepth[i],
				aerogelGreenIndex[i]);

		ss.str("");
		ss << "Aerogel" << i << "_Phys";
		new G4PVPlacement(0, G4ThreeVector(0 * cm, 0 * cm, aerogelZ[i]),
				ss.str().c_str(), aerogelLog, airPhys, false, 0);
	}

}

G4LogicalVolume* Cerenkov_Counter3::buildAerogelTile(const G4String& name,
		double length, double width, double depth, double greenIndex) {
	const double lambda0 = 83.22 * nm; // from Bellunato paper
	double a0 = scale_a0(greenIndex, lambda0); // a0 = 0.05639 -- from Bellunato paper
	const double scatLength_400 = 5.4 * cm;
	const int M = 21; // from quantum efficiency plot data points.

	double wavelength[M], aerogelRefIndex[M], aerogelEnergy[M],
			aerogelAbsorp[M], aerogelRayleighAbsorp[M];

	for (int i = 0; i < M; i++) {
		wavelength[i] = (200 * nm) + i * (25 * nm); // wavelengths range 275nm - 700nm
		aerogelEnergy[i] = twopi * hbarc / (wavelength[i]);
		aerogelRefIndex[i] = getIndex(wavelength[i], a0, lambda0);
		aerogelAbsorp[i] = scatLength_400
				* pow(((wavelength[i] / nm) / 400), 4);
		aerogelRayleighAbsorp[i] = 5 * cm;
	}

	G4Material* aerogel = new G4Material("Aerogel", 0.2 * g / cm3, 2);
	aerogel->AddMaterial(air, 0.35 * perCent);
	aerogel->AddMaterial(SiO2, 99.65 * perCent);

	G4MaterialPropertiesTable* aerogelMPT = new G4MaterialPropertiesTable();
	aerogelMPT->AddProperty("RINDEX", aerogelEnergy, aerogelRefIndex, M);
	aerogelMPT->AddProperty("ABSLENGTH", aerogelEnergy, aerogelAbsorp, M);
	aerogelMPT->AddProperty("RAYLEIGH", aerogelEnergy, aerogelRayleighAbsorp,
			M);
	aerogel->SetMaterialPropertiesTable(aerogelMPT);

	G4Box* aerogelSolid = new G4Box(
			name.substr(0, name.find('_', 1)) + "_solid", 0.5 * length,
			0.5 * width, 0.5 * depth);

	return new G4LogicalVolume(aerogelSolid, aerogel, name);
}

double Cerenkov_Counter3::scale_a0(double greenIndex, double lambda0) {
	return (pow(greenIndex, 2.) - 1) * (pow(532., 2.) - pow(lambda0 / nm, 2.))
			/ pow(532., 2.);
}

double Cerenkov_Counter3::getIndex(double lambda, double a0, double lambda0) {
	return sqrt(
			1
					+ a0 * pow(lambda / nm, 2.)
							/ (pow(lambda / nm, 2.) - pow(lambda0 / nm, 2.)));
}
